library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
#pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)
dataBreast <- read.csv("~/GitHub/RISKPLOTS/DATA/wpbc.data", header=FALSE)
table(dataBreast$V2)
##
## N R
## 151 47
rownames(dataBreast) <- dataBreast$V1
dataBreast$V1 <- NULL
dataBreast$status <- 1*(dataBreast$V2=="R")
dataBreast$V2 <- NULL
dataBreast$time <- dataBreast$V3
dataBreast$V3 <- NULL
dataBreast <- sapply(dataBreast,as.numeric)
## Warning in lapply(X = X, FUN = FUN, ...): NAs introduced by coercion
dataBreast <- as.data.frame(dataBreast[complete.cases(dataBreast),])
table(dataBreast$status)
##
## 0 1
## 148 46
ml <- BSWiMS.model(Surv(time,status)~1,data=dataBreast)
[++++++]
sm <- summary(ml)
pander::pander(sm$coefficients)
| Estimate | lower | HR | upper | u.Accuracy | r.Accuracy | |
|---|---|---|---|---|---|---|
| V24 | 4.69e-02 | 1.01 | 1.05 | 1.08 | 0.598 | 0.237 |
| V26 | 4.72e-03 | 1.00 | 1.00 | 1.01 | 0.593 | 0.237 |
| V27 | 2.42e-04 | 1.00 | 1.00 | 1.00 | 0.608 | 0.237 |
| V34 | 1.19e-02 | 1.00 | 1.01 | 1.02 | 0.634 | 0.237 |
| V7 | 6.05e-08 | 1.00 | 1.00 | 1.00 | 0.588 | 0.237 |
| V35 | 5.06e-06 | 1.00 | 1.00 | 1.00 | 0.727 | 0.237 |
| full.Accuracy | u.AUC | r.AUC | full.AUC | IDI | NRI | z.IDI | |
|---|---|---|---|---|---|---|---|
| V24 | 0.598 | 0.609 | 0.5 | 0.609 | 0.0619 | 0.437 | 2.87 |
| V26 | 0.593 | 0.598 | 0.5 | 0.598 | 0.0626 | 0.393 | 2.77 |
| V27 | 0.608 | 0.608 | 0.5 | 0.608 | 0.0563 | 0.434 | 2.76 |
| V34 | 0.634 | 0.618 | 0.5 | 0.618 | 0.0320 | 0.471 | 2.42 |
| V7 | 0.588 | 0.595 | 0.5 | 0.595 | 0.0487 | 0.380 | 2.30 |
| V35 | 0.727 | 0.641 | 0.5 | 0.641 | 0.0289 | 0.565 | 2.28 |
| z.NRI | Delta.AUC | Frequency | |
|---|---|---|---|
| V24 | 2.67 | 0.1091 | 1 |
| V26 | 2.38 | 0.0983 | 1 |
| V27 | 2.63 | 0.1084 | 1 |
| V34 | 2.85 | 0.1178 | 1 |
| V7 | 2.30 | 0.0949 | 1 |
| V35 | 3.50 | 0.1412 | 1 |
Here we evaluate the model using the RRPlot() function.
Here we will use the predicted event probability assuming a baseline hazard for events withing 5 years
index <- predict(ml,dataBreast)
timeinterval <- 2*mean(subset(dataBreast,status==1)$time)
h0 <- sum(dataBreast$status & dataBreast$time <= timeinterval)
h0 <- h0/sum((dataBreast$time > timeinterval) | (dataBreast$status==1))
pander::pander(t(c(h0=h0,timeinterval=timeinterval)),caption="Initial Parameters")
| h0 | timeinterval |
|---|---|
| 0.323 | 51.1 |
rdata <- cbind(dataBreast$status,ppoisGzero(index,h0))
rownames(rdata) <- rownames(dataBreast)
rrAnalysisTrain <- RRPlot(rdata,atRate=c(0.90),
timetoEvent=dataBreast$time,
title="Raw Train: Breast Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)
As we can see the Observed probability as well as the Time vs. Events are not calibrated.
pander::pander(t(rrAnalysisTrain$keyPoints),caption="Threshold values")
| @:0.9 | @MAX_BACC | @MAX_RR | @SPE100 | p(0.5) | |
|---|---|---|---|---|---|
| Thr | 0.42042 | 0.2640 | 0.1655 | 1.61e-01 | 0.497448 |
| RR | 2.18301 | 2.2857 | 4.3220 | 2.50e+01 | 2.307692 |
| RR_LCI | 1.30105 | 1.3037 | 0.6348 | 5.22e-02 | 1.275480 |
| RR_UCI | 3.66282 | 4.0074 | 29.4258 | 1.20e+04 | 4.175246 |
| SEN | 0.26087 | 0.6957 | 0.9783 | 1.00e+00 | 0.152174 |
| SPE | 0.89865 | 0.5608 | 0.1081 | 6.76e-02 | 0.952703 |
| BACC | 0.57976 | 0.6282 | 0.5432 | 5.34e-01 | 0.552438 |
| NetBenefit | 0.00578 | 0.0448 | 0.0971 | 1.00e-01 | 0.000374 |
pander::pander(t(rrAnalysisTrain$OERatio$estimate),caption="O/E Ratio")
| O/E | Low | Upper | p.value |
|---|---|---|---|
| 0.844 | 0.618 | 1.13 | 0.278 |
pander::pander(t(rrAnalysisTrain$OE95ci),caption="O/E Mean")
| mean | 50% | 2.5% | 97.5% |
|---|---|---|---|
| 1.02 | 1.02 | 0.976 | 1.08 |
pander::pander(t(rrAnalysisTrain$OAcum95ci),caption="O/Acum Mean")
| mean | 50% | 2.5% | 97.5% |
|---|---|---|---|
| 0.792 | 0.793 | 0.786 | 0.799 |
pander::pander(rrAnalysisTrain$c.index$cstatCI,caption="C. Index")
| mean.C Index | median | lower | upper |
|---|---|---|---|
| 0.68 | 0.678 | 0.594 | 0.754 |
pander::pander(t(rrAnalysisTrain$ROCAnalysis$aucs),caption="ROC AUC")
| est | lower | upper |
|---|---|---|
| 0.637 | 0.546 | 0.728 |
pander::pander((rrAnalysisTrain$ROCAnalysis$sensitivity),caption="Sensitivity")
| est | lower | upper |
|---|---|---|
| 0.261 | 0.143 | 0.411 |
pander::pander((rrAnalysisTrain$ROCAnalysis$specificity),caption="Specificity")
| est | lower | upper |
|---|---|---|
| 0.899 | 0.838 | 0.942 |
pander::pander(t(rrAnalysisTrain$thr_atP),caption="Probability Thresholds")
| 90% |
|---|
| 0.418 |
pander::pander(rrAnalysisTrain$surdif,caption="Logrank test")
| N | Observed | Expected | (O-E)^2/E | (O-E)^2/V | |
|---|---|---|---|---|---|
| class=0 | 167 | 34 | 41.1 | 1.23 | 11.6 |
| class=1 | 27 | 12 | 4.9 | 10.27 | 11.6 |
op <- par(no.readonly = TRUE)
calprob <- CoxRiskCalibration(ml,dataBreast,"status","time")
pander::pander(c(h0=calprob$h0,
Gain=calprob$hazardGain,
DeltaTime=calprob$timeInterval),
caption="Cox Calibration Parameters")
| h0 | Gain | DeltaTime |
|---|---|---|
| 0.25 | 0.774 | 40.9 |
h0 <- calprob$h0
timeinterval <- calprob$timeInterval;
rdata <- cbind(dataBreast$status,calprob$prob)
rrAnalysisTrain <- RRPlot(rdata,atRate=c(0.90),
timetoEvent=dataBreast$time,
title="Calibrated Train: Breast",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)
pander::pander(t(rrAnalysisTrain$keyPoints),caption="Threshold values")
| @:0.9 | @MAX_BACC | @MAX_RR | @SPE100 | p(0.5) | |
|---|---|---|---|---|---|
| Thr | 0.3445 | 0.2113 | 0.131 | 1.27e-01 | 0.49512 |
| RR | 2.1830 | 2.2857 | 4.322 | 2.50e+01 | 2.54422 |
| RR_LCI | 1.3011 | 1.3037 | 0.635 | 5.22e-02 | 1.27021 |
| RR_UCI | 3.6628 | 4.0074 | 29.426 | 1.20e+04 | 5.09603 |
| SEN | 0.2609 | 0.6957 | 0.978 | 1.00e+00 | 0.08696 |
| SPE | 0.8986 | 0.5608 | 0.108 | 6.76e-02 | 0.97973 |
| BACC | 0.5798 | 0.6282 | 0.543 | 5.34e-01 | 0.53334 |
| NetBenefit | 0.0212 | 0.0752 | 0.130 | 1.33e-01 | 0.00546 |
pander::pander(t(rrAnalysisTrain$OERatio$estimate),caption="O/E Ratio")
| O/E | Low | Upper | p.value |
|---|---|---|---|
| 0.87 | 0.637 | 1.16 | 0.408 |
pander::pander(t(rrAnalysisTrain$OE95ci),caption="O/E Mean")
| mean | 50% | 2.5% | 97.5% |
|---|---|---|---|
| 1.06 | 1.06 | 1.01 | 1.11 |
pander::pander(t(rrAnalysisTrain$OAcum95ci),caption="O/Acum Mean")
| mean | 50% | 2.5% | 97.5% |
|---|---|---|---|
| 0.962 | 0.962 | 0.953 | 0.97 |
pander::pander(t(rrAnalysisTrain$c.index$cstatCI),caption="C. Index")
| mean.C Index | median | lower | upper |
|---|---|---|---|
| 0.68 | 0.682 | 0.603 | 0.758 |
pander::pander(t(rrAnalysisTrain$ROCAnalysis$aucs),caption="ROC AUC")
| est | lower | upper |
|---|---|---|
| 0.637 | 0.546 | 0.728 |
pander::pander((rrAnalysisTrain$ROCAnalysis$sensitivity),caption="Sensitivity")
| est | lower | upper |
|---|---|---|
| 0.261 | 0.143 | 0.411 |
pander::pander((rrAnalysisTrain$ROCAnalysis$specificity),caption="Specificity")
| est | lower | upper |
|---|---|---|
| 0.899 | 0.838 | 0.942 |
pander::pander(t(rrAnalysisTrain$thr_atP),caption="Probability Thresholds")
| 90% |
|---|
| 0.343 |
pander::pander(rrAnalysisTrain$surdif,caption="Logrank test")
| N | Observed | Expected | (O-E)^2/E | (O-E)^2/V | |
|---|---|---|---|---|---|
| class=0 | 167 | 34 | 41.1 | 1.23 | 11.6 |
| class=1 | 27 | 12 | 4.9 | 10.27 | 11.6 |
Here we use the estimated h0 and timeinterval from the full set
rcv <- randomCV(theData=dataBreast,
theOutcome = Surv(time,status)~1,
fittingFunction=BSWiMS.model,
trainFraction = 0.9,
repetitions=100,
classSamplingType = "Pro"
)
.[+++++].[++++].[+++].[++].[++++++].[+++++++++].[-++].[+++].[+++].[++++]10 Tested: 127 Avg. Selected: 4.6 Min Tests: 1 Max Tests: 5 Mean Tests: 1.574803 . MAD: 0.4847067
.[++++++].[++++].[++++++++].[++++++++].[++++++++].[+].[+++++++].[+++++++].[–].[++++]20 Tested: 172 Avg. Selected: 5.3 Min Tests: 1 Max Tests: 6 Mean Tests: 2.325581 . MAD: 0.4774552
.[+++].[+++++].[+++].[+++].[++].[++++].[++++].[+++++].[-+++].[+++++]30 Tested: 187 Avg. Selected: 4.9 Min Tests: 1 Max Tests: 8 Mean Tests: 3.208556 . MAD: 0.4783071
.[++].[+++++].[+].[+++++].[+++].[++++].[++++++].[+++++].[+++].[++++]40 Tested: 191 Avg. Selected: 4.875 Min Tests: 1 Max Tests: 10 Mean Tests: 4.188482 . MAD: 0.4808877
.[++++++].[+++].[++++].[++++].[++++].[+++++++].[++++].[+++++].[+++++].[++++]50 Tested: 193 Avg. Selected: 4.98 Min Tests: 1 Max Tests: 12 Mean Tests: 5.181347 . MAD: 0.4835019
.[++++].[++++++].[++++].[+++++].[+++++].[++++++].[+++++].[++++].[+++].[+++]60 Tested: 194 Avg. Selected: 5.066667 Min Tests: 1 Max Tests: 13 Mean Tests: 6.185567 . MAD: 0.4835778
.[+++++++].[+++++++++].[++++++].[++++++].[+++++++].[+].[++++++++].[+++].[++++].[++++]70 Tested: 194 Avg. Selected: 5.2 Min Tests: 1 Max Tests: 16 Mean Tests: 7.216495 . MAD: 0.484655
.[++++].[++++].[++++].[+++].[+++].[+++].[++].[+++++++].[+++++].[+++++++]80 Tested: 194 Avg. Selected: 5.1625 Min Tests: 2 Max Tests: 16 Mean Tests: 8.247423 . MAD: 0.4841967
.[+++++].[+++++].[+++++].[++].[+++++].[+++++++].[++++].[+++].[+++++++].[++++]90 Tested: 194 Avg. Selected: 5.144444 Min Tests: 2 Max Tests: 17 Mean Tests: 9.278351 . MAD: 0.4839212
.[+++].[+++].[++++].[+++++].[+++++].[++++].[++++++].[+++].[+++++++].[+++]100 Tested: 194 Avg. Selected: 5.16 Min Tests: 2 Max Tests: 18 Mean Tests: 10.30928 . MAD: 0.4843842
stp <- rcv$survTestPredictions
stp <- stp[!is.na(stp[,4]),]
bbx <- boxplot(unlist(stp[,1])~rownames(stp),plot=FALSE)
times <- bbx$stats[3,]
status <- boxplot(unlist(stp[,2])~rownames(stp),plot=FALSE)$stats[3,]
prob <- ppoisGzero(boxplot(unlist(stp[,4])~rownames(stp),plot=FALSE)$stats[3,],h0)
rdatacv <- cbind(status,prob)
rownames(rdatacv) <- bbx$names
names(times) <- bbx$names
rrAnalysisTest <- RRPlot(rdatacv,atThr = rrAnalysisTrain$thr_atP,
timetoEvent=times,
title="Test: Breast Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)
pander::pander(t(rrAnalysisTest$keyPoints),caption="Threshold values")
| @:0.343 | @MAX_BACC | @MAX_RR | @SPE100 | p(0.5) | |
|---|---|---|---|---|---|
| Thr | 0.34323 | 0.2008 | 0.131 | 1.13e-01 | 0.49302 |
| RR | 1.40625 | 2.1967 | 2.814 | 1.47e+01 | 1.95767 |
| RR_LCI | 0.78018 | 1.2343 | 0.733 | 3.13e-02 | 0.89980 |
| RR_UCI | 2.53472 | 3.9096 | 10.809 | 6.89e+03 | 4.25926 |
| SEN | 0.21739 | 0.7174 | 0.957 | 1.00e+00 | 0.08696 |
| SPE | 0.85135 | 0.5203 | 0.135 | 4.05e-02 | 0.96622 |
| BACC | 0.53437 | 0.6188 | 0.546 | 5.20e-01 | 0.52659 |
| NetBenefit | -0.00771 | 0.0781 | 0.128 | 1.44e-01 | -0.00444 |
pander::pander(t(rrAnalysisTest$OERatio$estimate),caption="O/E Ratio")
| O/E | Low | Upper | p.value |
|---|---|---|---|
| 0.849 | 0.622 | 1.13 | 0.307 |
pander::pander(t(rrAnalysisTest$OE95ci),caption="O/E Mean")
| mean | 50% | 2.5% | 97.5% |
|---|---|---|---|
| 1.02 | 1.02 | 0.973 | 1.08 |
pander::pander(t(rrAnalysisTest$OAcum95ci),caption="O/Acum Mean")
| mean | 50% | 2.5% | 97.5% |
|---|---|---|---|
| 0.869 | 0.869 | 0.852 | 0.885 |
pander::pander(rrAnalysisTest$c.index$cstatCI,caption="C. Index")
| mean.C Index | median | lower | upper |
|---|---|---|---|
| 0.649 | 0.647 | 0.565 | 0.728 |
pander::pander(t(rrAnalysisTest$ROCAnalysis$aucs),caption="ROC AUC")
| est | lower | upper |
|---|---|---|
| 0.603 | 0.509 | 0.696 |
pander::pander((rrAnalysisTest$ROCAnalysis$sensitivity),caption="Sensitivity")
| est | lower | upper |
|---|---|---|
| 0.217 | 0.109 | 0.364 |
pander::pander((rrAnalysisTest$ROCAnalysis$specificity),caption="Specificity")
| est | lower | upper |
|---|---|---|
| 0.851 | 0.784 | 0.904 |
pander::pander(t(rrAnalysisTest$thr_atP),caption="Probability Thresholds")
| 90% |
|---|
| 0.343 |
pander::pander(rrAnalysisTest$surdif,caption="Logrank test")
| N | Observed | Expected | (O-E)^2/E | (O-E)^2/V | |
|---|---|---|---|---|---|
| class=0 | 162 | 36 | 39.73 | 0.351 | 2.6 |
| class=1 | 32 | 10 | 6.27 | 2.225 | 2.6 |
rdatacv <- cbind(status,prob,times)
calprob <- CalibrationProbPoissonRisk(rdatacv)
pander::pander(c(h0=calprob$h0,
Gain=calprob$hazardGain,
DeltaTime=calprob$timeInterval),
caption="Cox Calibration Parameters")
| h0 | Gain | DeltaTime |
|---|---|---|
| 0.316 | 0.979 | 41.4 |
timeinterval <- calprob$timeInterval;
rdata <- cbind(status,calprob$prob)
rrAnalysisTest <- RRPlot(rdata,atRate=c(0.90),
timetoEvent=times,
title="Calibrated Test: Breast",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)
### Calibrated Test Performance
pander::pander(t(rrAnalysisTest$keyPoints),caption="Threshold values")
| @:0.9 | @MAX_BACC | @MAX_RR | @SPE100 | p(0.5) | |
|---|---|---|---|---|---|
| Thr | 0.3620 | 0.1970 | 0.128 | 1.11e-01 | 5.01e-01 |
| RR | 1.3344 | 2.1967 | 2.814 | 1.47e+01 | 2.21e+00 |
| RR_LCI | 0.6783 | 1.2343 | 0.733 | 3.13e-02 | 1.05e+00 |
| RR_UCI | 2.6251 | 3.9096 | 10.809 | 6.89e+03 | 4.65e+00 |
| SEN | 0.1522 | 0.7174 | 0.957 | 1.00e+00 | 8.70e-02 |
| SPE | 0.8919 | 0.5203 | 0.135 | 4.05e-02 | 9.73e-01 |
| BACC | 0.5220 | 0.6188 | 0.546 | 5.20e-01 | 5.30e-01 |
| NetBenefit | -0.0107 | 0.0803 | 0.130 | 1.46e-01 | -8.43e-05 |
pander::pander(t(rrAnalysisTest$OERatio$estimate),caption="O/E Ratio")
| O/E | Low | Upper | p.value |
|---|---|---|---|
| 0.876 | 0.641 | 1.17 | 0.407 |
pander::pander(t(rrAnalysisTest$OE95ci),caption="O/E Mean")
| mean | 50% | 2.5% | 97.5% |
|---|---|---|---|
| 1.06 | 1.06 | 1 | 1.11 |
pander::pander(t(rrAnalysisTest$OAcum95ci),caption="O/Acum Mean")
| mean | 50% | 2.5% | 97.5% |
|---|---|---|---|
| 0.883 | 0.884 | 0.867 | 0.9 |
pander::pander(rrAnalysisTest$c.index$cstatCI,caption="C. Index")
| mean.C Index | median | lower | upper |
|---|---|---|---|
| 0.649 | 0.648 | 0.568 | 0.723 |
pander::pander(t(rrAnalysisTest$ROCAnalysis$aucs),caption="ROC AUC")
| est | lower | upper |
|---|---|---|
| 0.603 | 0.509 | 0.696 |
pander::pander((rrAnalysisTest$ROCAnalysis$sensitivity),caption="Sensitivity")
| est | lower | upper |
|---|---|---|
| 0.152 | 0.0634 | 0.289 |
pander::pander((rrAnalysisTest$ROCAnalysis$specificity),caption="Specificity")
| est | lower | upper |
|---|---|---|
| 0.899 | 0.838 | 0.942 |
pander::pander(t(rrAnalysisTest$thr_atP),caption="Probability Thresholds")
| 90% |
|---|
| 0.364 |
pander::pander(rrAnalysisTest$surdif,caption="Logrank test")
| N | Observed | Expected | (O-E)^2/E | (O-E)^2/V | |
|---|---|---|---|---|---|
| class=0 | 172 | 39 | 42.02 | 0.217 | 2.54 |
| class=1 | 22 | 7 | 3.98 | 2.295 | 2.54 |